In this script we created different plots for the description of the NBI population. We created Timelines, Population size per year and Barplots for all individuals and the alive females to show the distribution between the sexes, stages, colonies, raising types and the mortality types.

Setup

## libraries used in this script
library("here")
library("tidyverse")
library("viridisLite")

Data

NBI <- readRDS(file = here::here("output", "data-proc", "NBI_data.rds"))

Save some standard variables

mar_default <- c(5, 4, 4, 2) + 0.1
mar_small <- c(0,0,0,0)
zero_line <- rep(0, 12)

Population size per year

Plot the whole Population size per year

# Calculate the number of living individuals per year
popsize_NBI <- tibble(
  "2008" = ifelse(test = NBI$hatch_year <= 2008 & (NBI$death_year >= 2008 | NBI$death_year == 0), 
                 yes = 1, 
                 no = 0) %>% 
          sum(), 
  "2009" = ifelse(test = NBI$hatch_year <= 2009 & (NBI$death_year >= 2009 | NBI$death_year == 0), 
                 yes = 1, 
                 no = 0) %>% 
          sum(),
  "2010" = ifelse(test = NBI$hatch_year <= 2010 & (NBI$death_year >= 2010 | NBI$death_year == 0), 
                 yes = 1, 
                 no = 0) %>% 
          sum(),
  "2011" = ifelse(test = NBI$hatch_year <= 2011 & (NBI$death_year >= 2011 | NBI$death_year == 0), 
                 yes = 1, 
                 no = 0) %>% 
          sum(),
  "2012" = ifelse(test = NBI$hatch_year <= 2012 & (NBI$death_year >= 2012 | NBI$death_year == 0), 
                 yes = 1, 
                 no = 0) %>% 
          sum(),
  "2013" = ifelse(test = NBI$hatch_year <= 2013 & (NBI$death_year >= 2013 | NBI$death_year == 0), 
                 yes = 1, 
                 no = 0) %>% 
          sum(),
  "2014" = ifelse(test = NBI$hatch_year <= 2014 & (NBI$death_year >= 2014 | NBI$death_year == 0), 
                 yes = 1, 
                 no = 0) %>% 
          sum(),
  "2015" = ifelse(test = NBI$hatch_year <= 2015 & (NBI$death_year >= 2015 | NBI$death_year == 0), 
                 yes = 1, 
                 no = 0) %>% 
          sum(),
  "2016" = ifelse(test = NBI$hatch_year <= 2016 & (NBI$death_year >= 2016 | NBI$death_year == 0), 
                 yes = 1, 
                 no = 0) %>% 
          sum(),
  "2017" = ifelse(test = NBI$hatch_year <= 2017 & (NBI$death_year >= 2017 | NBI$death_year == 0), 
                 yes = 1, 
                 no = 0) %>% 
          sum(),
  "2018" = ifelse(test = NBI$hatch_year <= 2018 & (NBI$death_year >= 2018 | NBI$death_year == 0), 
                 yes = 1, 
                 no = 0) %>% 
          sum(),
  "2019" = ifelse(test = NBI$hatch_year <= 2019 & (NBI$death_year >= 2019 | NBI$death_year == 0), 
                 yes = 1, 
                 no = 0) %>% 
          sum())

# Calculate the change per year
for (i in 2:ncol(popsize_NBI)) {
  
  popsize_NBI[2,i] <- popsize_NBI[1,i] - popsize_NBI[1,i-1]

}
# only numbers
pop_num <- as.numeric(popsize_NBI[1,])

# plot it
par(mar = (c(5, 4.5, 1, 2) + 0.1))

plot(pop_num, type = "n", xaxt = "n", yaxt = "n", xlab = "Year", ylab = "Number of individuals", cex.lab = 2, ylim = c(0,200))
axis(1, at = 1:12, labels = 2008:2019, cex.axis = 2)
axis(2, cex.axis = 2)

lines(zero_line, col = "grey85")
polygon(c(1:12, rev(1:12)), c(pop_num, rev(zero_line)), col = "grey85", border= NA)

lines(pop_num, lwd = 3)
points(pop_num, pch = 16, cex = 1.5)

box()

par(mar = mar_default)

Save the plot

# png(filename = here("plots", "01_description", "Popsize_per_year.png"))
# par(mar = (c(5, 4.5, 1, 2) + 0.1))
# 
# plot(pop_num, type = "n", xaxt = "n", yaxt = "n", xlab = "Year", ylab = "Number of individuals", cex.lab = 2, ylim = c(0,200))
# axis(1, at = 1:12, labels = 2008:2019, cex.axis = 2)
# axis(2, cex.axis = 2)
# 
# lines(zero_line, col = "grey85")
# polygon(c(1:12, rev(1:12)), c(pop_num, rev(zero_line)), col = "grey85", border= NA)
# 
# lines(pop_num, lwd = 3)
# points(pop_num, pch = 16, cex = 1.5)
# 
# box()
# par(mar = mar_default)
# 
# dev.off()

Customized barplots

Overall population

Show the values for each category and build a matrix with them. The NAs between some values are there for some space between the categories

# Comparison between delivered individuals and all other individuals
NBI_deliv <- NBI[NBI$mortality_cause == "delivered",] # 9 delivered individuals (delivered to a colony outside the project)
NBI_proj <- NBI[NBI$mortality_cause != "delivered",] # 375 not delivered individuals

# Comparison of the reached stages and the life status (alive or dead)
table(NBI_proj$life_status, NBI_proj$reached_stage)
##    
##       1   2   3   4
##   0  71  28  13  31
##   1 123  45  33  31
table(NBI_deliv$life_status, NBI_deliv$reached_stage) # all alive at the time they were delivered 
##    
##     1 2 3 4
##   0 8 1 0 0
##   1 0 0 0 0
# List the number of individuals for different characteristics
table(NBI$reached_stage)
## 
##   1   2   3   4 
## 202  74  46  62
table(NBI$sex)
## 
##   f   m   u 
## 184 195   5
table(NBI$breeding_area)
## 
##  Burghausen       Kuchl  allocation Ueberlingen      Rosegg 
##         127          89          57          89          22
table(NBI$raising_type)
## 
##   foster parent          parent supplementation 
##             213             162               9
table(NBI$mortality_cause)
## 
##             alive         delivered during human care     electrocution    injury/illness             other         predation              shot           unknown 
##               143                 9                 4                46                23                 2                15                18               124
NBI_mat <- matrix(data = c(
             202, 74,46,62,NA,           # Stage
             184,195, 5,NA,              # Sex
             127, 89,89,22,57,NA,        # Colony (B, K, Ue, R, Allo)
             213, 162,9, NA,             # Raising type
             4,46,23,15,18,2,124,9,143)) # Mortality  type (d h c, elec, inj, pred, shot. other, unk, del, alive)

Plot

par(mar = c(5,5,2,0))
barplot(NBI_mat, beside = T, ylim = c(0,220), ylab = "Number of Individuals", border = NA, 
        cex.axis = 3, cex.lab = 3,
        col = c(
                magma(4, begin = 0.25, end = 0.6), "white", # Stage
                
                plasma(1, direction = -1, begin = 0.0, end = 0.6), viridis(1, begin = 0.2, alpha = 0.9), # Sex
                  magma(1, begin = 0.1, alpha = 0.6),"white", # Sex
                
                inferno(5,begin =0.65,  end = 0.95, direction = -1), "white", # Colony 
                
                viridis(3, begin = 0.65, end = 0.85),"white", # Raising type
                
                gray(0:7 / 8), viridis(1, begin = 0.5))) # Mortality type

axis(1, at = c(3,7.5,12.5,17.5,24.5), labels = c("", "", "", "", ""), cex.axis = 3)
axis(1, at = c(3,7.5,12.5,17.5,24.5), labels = c("Stage", "Sex", "Colony", "Raising type", "Mortality type"), cex.axis = 3, line = 1, lwd = 0)
box()

Save the Plot

# png(filename = here("plots", "01_description", "Barplot_all.png"), width = 1300, height = 835)
# par(mar = c(5,5,2,0))
# 
# barplot(NBI_mat, beside = T, ylim = c(0,220), ylab = "Number of Individuals", border = NA, 
#         cex.axis = 3, cex.lab = 3,
#         col = c(
#                 magma(4, begin = 0.25, end = 0.6), "white", # Stage
#                 
#                 plasma(1, direction = -1, begin = 0.0, end = 0.6), viridis(1, begin = 0.2, alpha = 0.9), # Sex
#                   magma(1, begin = 0.1, alpha = 0.6),"white", # Sex
#                 
#                 inferno(5,begin =0.65,  end = 0.95, direction = -1), "white", # Colony 
#                 
#                 viridis(3, begin = 0.65, end = 0.85),"white", # Raising type
#                 
#                 gray(0:7 / 8), viridis(1, begin = 0.5))) # Mortality type
# 
# axis(1, at = c(3,7.5,12.5,17.5,24.5), labels = c("", "", "", "", ""), cex.axis = 3)
# axis(1, at = c(3,7.5,12.5,17.5,24.5), labels = c("Stage", "Sex", "Colony", "Raising type", "Mortality type"), cex.axis = 3, line = 1, lwd = 0)
# box()
# 
# dev.off()

Legend

par(mar = mar_small)
plot(1, type = "n", axes = "F", xlab = "", ylab = "") # empty plot

legend(x = "center", 
       legend = c("Stage","Juveniles","1 Year Old","2 Year Old","Adults",   "",
                  "Sex","Females","Males","Unknown",   "", 
                  "Colony","Burghausen","Kuchl","Überlingen", "Rosegg", "unassigned yet",   "", 
                  "Raising type","Foster parents","Biological Parents", "Supplementation",   "", 
                  "Mortality type","During human care","Electrocution","Injury or illness","Predation","Shot","Other", "Unknown", "Delivered", "Alive"), 
       fill = c("white",magma(4, begin = 0.25, end = 0.6),   "white", # Stage
                "white",plasma(1, direction = -1, begin = 0.0, end = 0.6), viridis(1, begin = 0.2, alpha = 0.9), # Sex
                        magma(1, begin = 0.1, alpha = 0.6),"white", # Sex
                "white",inferno(5,begin =0.65,  end = 0.95, direction = -1),   "white", # Colony
                "white",viridis(3, begin = 0.65, end = 0.85),         "white", # Raising type
                "white" , gray(0:7 / 8), viridis(1, begin = 0.5)), # Mortality type
       border = NA, cex = 3)

par(mar = mar_default)

Save the Legend

# png(filename = here("plots", "01_description", "Barplot_all_legend.png"), width = 2174, height = 1800)
# 
# par(mar = mar_small)
# plot(1, type = "n", axes = "F", xlab = "", ylab = "") # empty plot
# 
# legend(x = "center", 
#        legend = c("Stage","Juveniles","1 Year Old","2 Year Old","Adults",   "",
#                   "Sex","Females","Males","Unknown",   "", 
#                   "Colony","Burghausen","Kuchl","Überlingen", "Rosegg", "unassigned yet",   "", 
#                   "Raising type","Foster parents","Biological Parents", "Supplementation",   "", 
#                   "Mortality type","During human care","Electrocution","Injury or illness","Predation","Shot","Other", "Unknown", "Delivered", "Alive"), 
#        fill = c("white",magma(4, begin = 0.25, end = 0.6),   "white", # Stage
#                 "white",plasma(1, direction = -1, begin = 0.0, end = 0.6), viridis(1, begin = 0.2, alpha = 0.9), # Sex
#                         magma(1, begin = 0.1, alpha = 0.6),"white", # Sex
#                 "white",inferno(5,begin =0.65,  end = 0.95, direction = -1),   "white", # Colony
#                 "white",viridis(3, begin = 0.65, end = 0.85),         "white", # Raising type
#                 "white" , gray(0:7 / 8), viridis(1, begin = 0.5)), # Mortality type
#        border = NA, cex = 3)
# 
# par(mar = mar_default)
# 
# dev.off()

Alive female individuals

table(NBI$sex, NBI$project_status) # 74 alive females, 69 alive males = 143 alive individuals
##    
##       0   1
##   f  74 110
##   m  69 126
##   u   0   5
# subset the data
NBI_f <- subset(x = NBI, subset = NBI$sex == "f" & NBI$project_status == "0")
table(NBI_f$reached_stage)
## 
##  1  2  3  4 
## 37 11  8 18
table(NBI_f$breeding_area)
## 
##  Burghausen       Kuchl  allocation Ueberlingen      Rosegg 
##          19          15           6          28           6
table(NBI_f$raising_type)
## 
##   foster parent          parent supplementation 
##              41              32               1
NBI_mat_f <- matrix(data = c(
             37, 11,8,18,NA,           # Stage
             19, 15,28,6,6,NA,         # Colony (B, K, Ue, R, Allo)
             41, 32,1))                # Raising type

Plot

par(mar = c(5,5,2,0))
barplot(NBI_mat_f, beside = T, ylim = c(0,50), ylab = "Number of Individuals", border = NA, 
        cex.axis = 3, cex.lab = 3,
        col = c(
                magma(4, begin = 0.25, end = 0.6), "white", # Stage
                
                inferno(5,begin =0.65,  end = 0.95, direction = -1), "white", # Colony 
                
                viridis(3, begin = 0.65, end = 0.85) # Raising type
                
                ))
axis(1, at = c(3,8.5,13.5), labels = c("", "", ""), cex.axis = 3)
axis(1, at = c(3,8.5,13.5), labels = c("Stage", "Colony", "Raising type"), cex.axis = 3, line = 1, lwd = 0)

box()

Save the Plot

# png(filename = here("plots", "01_description", "Barplot_alive_f.png"), width = 1300, height = 835)
# par(mar = c(5,5,2,0))
# 
# barplot(NBI_mat_f, beside = T, ylim = c(0,50), ylab = "Number of Individuals", border = NA, 
#         cex.axis = 3, cex.lab = 3,
#         col = c(
#                 magma(4, begin = 0.25, end = 0.6), "white", # Stage
#                 
#                 inferno(5,begin =0.65,  end = 0.95, direction = -1), "white", # Colony 
#                 
#                 viridis(3, begin = 0.65, end = 0.85) # Raising type
#                 
#                 ))
# axis(1, at = c(3,8.5,13.5), labels = c("", "", ""), cex.axis = 3)
# axis(1, at = c(3,8.5,13.5), labels = c("Stage", "Colony", "Raising type"), cex.axis = 3, line = 1, lwd = 0)
# 
# box()
# 
# dev.off()

Legend

par(mar = mar_small)
plot(1, type = "n", axes = "F", xlab = "", ylab = "")
legend("center", legend = c("Stage","Juveniles","1 Year Old","2 Year Old","Adults",   "",
                            "Colony","Burghausen","Kuchl","Überlingen", "Rosegg", "unassigned yet",   "",
                            "Raising type","Foster parents","Biological Parents", "Supplementation"),
       fill = c("white",magma(4, begin = 0.25, end = 0.6),   "white",
                "white",inferno(5,begin =0.65,  end = 0.95, direction = -1),   "white",
                "white",viridis(3, begin = 0.65, end = 0.85)), border = NA, cex = 1)

par(mar = mar_default)

Save the legend

# png(filename = here("plots", "01_description", "Barplot_alive_f_legend.png"), res = 130) # width = 300, height = 400
# par(mar = mar_small)
# 
# plot(1, type = "n", axes = "F", xlab = "", ylab = "")
# legend("center", legend = c("Stage","Juveniles","1 Year Old","2 Year Old","Adults",   "",
#                             "Colony","Burghausen","Kuchl","Überlingen", "Rosegg", "unassigned yet",   "",
#                             "Raising type","Foster parents","Biological Parents", "Supplementation"),
#        fill = c("white",magma(4, begin = 0.25, end = 0.6),   "white",
#                 "white",inferno(5,begin =0.65,  end = 0.95, direction = -1),   "white",
#                 "white",viridis(3, begin = 0.65, end = 0.85)), border = NA, cex = 1)
# 
# par(mar = mar_default)
# dev.off()

Session Info

## DO NOT REMOVE!
## We store the settings of your computer and the current versions of the
## packages used to allow for reproducibility
Sys.time()
## [1] "2022-01-06 19:33:39 CET"
#git2r::repository() ## uncomment if you are using GitHub
sessionInfo()
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 22000)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252    LC_MONETARY=German_Germany.1252 LC_NUMERIC=C                    LC_TIME=German_Germany.1252    
## 
## attached base packages:
## [1] splines   stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] viridisLite_0.3.0  eeptools_1.2.4     forcats_0.5.0      stringr_1.4.0      dplyr_1.0.2        purrr_0.3.4        readr_1.4.0        tidyr_1.1.2        tibble_3.0.3      
## [10] ggplot2_3.3.2      tidyverse_1.3.0    pwr_1.3-0          plot3D_1.4         MuMIn_1.43.17      modelsummary_0.6.5 mgcv_1.8-33        nlme_3.1-149       lattice_0.20-41   
## [19] here_0.1           ggeffects_1.0.0    gam_1.20           foreach_1.5.1      effects_4.2-0      carData_3.0-4     
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.3.1        backports_1.1.10    Hmisc_4.4-1         sp_1.4-4            usethis_1.6.3       digest_0.6.25       htmltools_0.5.0     fansi_0.4.1        
##   [9] magrittr_1.5        checkmate_2.0.0     memoise_1.1.0       cluster_2.1.0       remotes_2.2.0       modelr_0.1.8        sysfonts_0.8.1      rmdformats_0.3.7   
##  [17] prettyunits_1.1.1   jpeg_0.1-8.1        colorspace_1.4-1    blob_1.2.1          rvest_0.3.6         mitools_2.4         haven_2.3.1         xfun_0.18          
##  [25] tcltk_4.0.3         callr_3.5.0         crayon_1.3.4        jsonlite_1.7.1      lme4_1.1-25         zoo_1.8-8           survival_3.2-7      iterators_1.0.13   
##  [33] glue_1.4.2          kableExtra_1.3.1    gtable_0.3.0        webshot_0.5.2       pkgbuild_1.1.0      abind_1.4-5         scales_1.1.1        DBI_1.1.0          
##  [41] Rcpp_1.0.5          showtextdb_3.0      htmlTable_2.1.0     foreign_0.8-80      Formula_1.2-4       stats4_4.0.3        survey_4.0          vcd_1.4-8          
##  [49] htmlwidgets_1.5.2   httr_1.4.2          RColorBrewer_1.1-2  ellipsis_0.3.1      pkgconfig_2.0.3     farver_2.0.3        nnet_7.3-14         dbplyr_1.4.4       
##  [57] utf8_1.1.4          tidyselect_1.1.0    labeling_0.3        rlang_0.4.7         munsell_0.5.0       cellranger_1.1.0    tools_4.0.3         cli_2.0.2          
##  [65] generics_0.0.2      sjlabelled_1.1.7    devtools_2.3.2      broom_0.7.4         evaluate_0.14       arm_1.11-2          yaml_2.2.1          tables_0.9.6       
##  [73] processx_3.4.4      knitr_1.30          fs_1.5.0            showtext_0.9        xml2_1.3.2          compiler_4.0.3      rstudioapi_0.11     png_0.1-7          
##  [81] testthat_2.3.2      reprex_0.3.0        statmod_1.4.35      stringi_1.5.3       highr_0.8           ps_1.4.0            desc_1.2.0          Matrix_1.2-18      
##  [89] nloptr_1.2.2.2      vctrs_0.3.4         pillar_1.4.6        lifecycle_0.2.0     lmtest_0.9-38       estimability_1.3    maptools_1.0-2      data.table_1.13.2  
##  [97] insight_0.11.1      R6_2.4.1            latticeExtra_0.6-29 bookdown_0.20       gridExtra_2.3       sessioninfo_1.1.1   codetools_0.2-16    boot_1.3-25        
## [105] MASS_7.3-53         assertthat_0.2.1    pkgload_1.1.0       rprojroot_1.3-2     withr_2.3.0         d6_0.1.0.0          hms_0.5.3           grid_4.0.3         
## [113] rpart_4.1-15        coda_0.19-4         minqa_1.2.4         rmarkdown_2.4       misc3d_0.9-0        git2r_0.27.1        lubridate_1.7.9     base64enc_0.1-3    
## [121] tinytex_0.26